#include <assert.h>
#include <stdio.h>
#include <stdlib.h>

#include "../src/distribution.c"
#include "../src/qfunc.c"

typedef struct _QFuncTestDatum QFuncTestDatum ;
struct _QFuncTestDatum 
{ 
  struct
    {
      float     kappa;
      float     logalpha;
      float     tau[(5-1)*4];
      float     loglambda[4];
      float     pz[5];
      Rating    ratings[8];
      float     logalphastddev;
    }                   inputs;
  struct
    {
      float q;
      float dqdlogalpha;
      float ddqddlogalpha;
    }                   desired_output;
};

static void
test_qfunc (void)
{
  QFuncTestDatum data[] = 
    {{{0.983886,  0.464123, {0.068732, 0.516637, 0.958455, 0.926766, 0.186644,  0.853853, 0.979339, 0.341647, 0.535515, 0.475623, 0.905311,  0.10016, 0.460451, 0.0443592, 0.757365, 0.303473}, {0.428895,  0.664885, 0.2102, 0.987802}, {0.481582, 0.154956, 0.226314,  0.523679,  0.41285}, {{3, 3}, {2, 1}, {0, 4}, {0, 1}, {3, 3}, {0, 0}, {3,  1}, {2, 0}},  0.638319}, {-96.1904, -86.1887, -95.9791}}, {{0.267859,  0.596913, {0.226206, 0.784466, 0.288521, 0.255266, 0.690691,  0.308843, 0.383209, 0.155106, 0.23024, 0.264484, 0.625844,  0.851633, 0.801345, 0.599598, 0.415644, 0.863831}, {0.319763,  0.444642, 0.189329, 0.340152}, {0.906913, 0.806323, 0.92147,  0.743239,  0.680707}, {{3, 0}, {3, 1}, {1, 1}, {2, 1}, {2, 4}, {0, 0}, {3,  1}, {0, 2}},  0.0218568}, {-258.232, -813.338, -2231.31}}, {{0.632949,  0.487974, {0.990016, 0.713014, 0.24974, 0.332868, 0.759776,  0.44853, 0.623896, 0.481235, 0.958431, 0.848932, 0.208252,  0.617404, 0.638668, 0.40429, 0.0189228, 0.277252}, {0.731755,  0.597967, 0.0974527, 0.534013}, {0.0510484, 0.57611, 0.464503,  0.0460396,  0.0610325}, {{0, 1}, {2, 1}, {2, 0}, {1, 1}, {1, 3}, {1, 3}, {1,  4}, {0, 1}},  0.863096}, {-45.0952, -37.7032, -43.0054}}, {{0.214763,  0.713172, {0.301257, 0.414566, 0.590867, 0.231937, 0.342826,  0.565634, 0.382614, 0.614533, 0.704157, 0.161344, 0.363692,  0.33728, 0.972402, 0.563377, 0.266239, 0.803267}, {0.921354,  0.987266, 0.801736, 0.757228}, {0.860321, 0.12417, 0.586973,  0.0440556,  0.559065}, {{3, 0}, {1, 2}, {2, 2}, {2, 3}, {0, 3}, {1, 3}, {0,  0}, {1, 0}},  0.709604}, {-190.215, -184.416, -189.597}}, {{0.996106,  0.812118, {0.216239, 0.14397, 0.613491, 0.197586, 0.512082,  0.982626, 0.2498, 0.860305, 0.53968, 0.419249, 0.983561,  0.0570381, 0.618326, 0.431983, 0.181825, 0.29981}, {0.758005,  0.307814, 0.594852, 0.255755}, {0.19894, 0.59821, 0.598746,  0.443636,  0.982701}, {{0, 2}, {3, 1}, {2, 4}, {1, 3}, {0, 4}, {1, 4}, {1,  4}, {3, 1}},  0.454241}, {-177.556, -163.168, -178.322}}, {{0.985255,  0.246051, {0.470619, 0.471615, 0.735455, 0.385745, 0.930939,  0.0523655, 0.751894, 0.328707, 0.312613, 0.620382, 0.570069,  0.0288968, 0.554609, 0.312568, 0.975217, 0.773142}, {0.355669,  0.714358, 0.37647, 0.329506}, {0.372968, 0.260118, 0.391215,  0.0834551,  0.902349}, {{0, 3}, {0, 3}, {1, 2}, {3, 3}, {1, 0}, {1, 4}, {1,  1}, {3, 1}},  0.788503}, {-82.7167, -69.0279, -78.0789}}, {{0.65576,  0.69771, {0.971409, 0.736137, 0.903866, 0.369003, 0.658796,  0.115755, 0.333796, 0.340106, 0.104187, 0.803187, 0.35858,  0.566964, 0.748518, 0.0888285, 0.982109, 0.237458}, {0.375551,  0.828711, 0.590894, 0.154003}, {0.473202, 0.0402083, 0.935134,  0.456293,  0.501793}, {{0, 4}, {3, 2}, {3, 3}, {2, 4}, {2, 4}, {1, 3}, {0,  1}, {2, 0}},  0.304071}, {-154.337, -143.53, -162.312}}, {{0.0312683,  0.0872907, {0.842997, 0.188316, 0.697472, 0.747185, 0.73881,  0.385129, 0.338892, 0.180221, 0.990291, 0.296301, 0.356783,  0.942763, 0.614741, 0.46759, 0.765889, 0.78876}, {0.141539,  0.427382, 0.830755, 0.332467}, {0.639746, 0.123311, 0.799486,  0.245176,  0.79675}, {{0, 2}, {2, 0}, {3, 0}, {0, 1}, {3, 2}, {0, 0}, {3,  4}, {3, 3}},  0.934995}, {-81.1802, -62.2368, -71.8235}}, {{0.102014,  0.497991, {0.05794, 0.549865, 0.763122, 0.31777, 0.0676487,  0.253564, 0.40634, 0.375007, 0.452908, 0.785974, 0.640451,  0.586246, 0.311369, 0.358593, 0.809696, 0.253779}, {0.671622,  0.235282, 0.0102098, 0.00860281}, {0.874873, 0.300287, 0.908195,  0.510611,  0.816933}, {{0, 0}, {1, 0}, {3, 4}, {3, 0}, {3, 3}, {2, 4}, {1,  4}, {2, 3}},  0.750422}, {-173.604, -151.065, -163.186}}, {{0.145073,  0.192841, {0.749284, 0.496858, 0.738733, 0.817835, 0.296376,  0.710884, 0.0982826, 0.231588, 0.985007, 0.352291, 0.288587,  0.977809, 0.313385, 0.117009, 0.278377, 0.969206}, {0.438512,  0.816722, 0.370181, 0.458595}, {0.621579, 0.0662993, 0.225109,  0.265753,  0.872295}, {{3, 3}, {0, 3}, {3, 2}, {3, 1}, {2, 4}, {1, 2}, {0,  3}, {3, 4}}, 0.569441}, {-83.0251, -71.8233, -82.0265}}};

  for (unsigned int i = 0; i < sizeof (data) / sizeof (data[0]); ++i)
    {
      ParameterizedNormalDistribution prior_logalpha = 
        parameterized_normal_distribution (data[i].inputs.logalphastddev);

      float hq = 0;

      for (unsigned int k = 0; k < 5; ++k)
        {
          hq -= data[i].inputs.pz[k] * log (data[i].inputs.pz[k]);
        }

      float dqdlogalpha;
      float ddqddlogalpha;

      float q = qfunc (data[i].inputs.kappa,
                       data[i].inputs.logalpha,
                       &dqdlogalpha,
                       &ddqddlogalpha,
                       5,
                       data[i].inputs.tau,
                       data[i].inputs.loglambda,
                       data[i].inputs.pz,
                       data[i].inputs.ratings,
                       8,
                       &prior_logalpha.base,
                       hq);

      assert (fabs (q - data[i].desired_output.q) <=
              1e-4 * (1 + fabs (q) + fabs (data[i].desired_output.q)) ||
              (fprintf (stderr, "%g ?= %g\n", q, data[i].desired_output.q),
               0));

      assert (fabs (dqdlogalpha - data[i].desired_output.dqdlogalpha) <=
              5e-5 * (1 + fabs (dqdlogalpha) + fabs (data[i].desired_output.dqdlogalpha)) ||
              (fprintf (stderr, "%g ?= %g\n", dqdlogalpha, data[i].desired_output.dqdlogalpha),
               0));

      assert (fabs (ddqddlogalpha - data[i].desired_output.ddqddlogalpha) <=
              5e-5 * (1 + fabs (ddqddlogalpha) + fabs (data[i].desired_output.ddqddlogalpha)) ||
              (fprintf (stderr, "%g ?= %g\n", ddqddlogalpha, data[i].desired_output.ddqddlogalpha),
               0));
    }
}

typedef struct _DQFuncTestDatum DQFuncTestDatum ;
struct _DQFuncTestDatum 
{ 
  struct
    {
      float     kappa;
      float     logalpha;
      float     tau[(5-1)*4];
      float     loglambda[4];
      float     pz[5];
      Rating    ratings[8];
      float     logalphastddev;
    }                   inputs;
  struct
    {
      float dtau[(5-1)*4];
      float dloglambda[4];
    }                   desired_output;
};

static void
test_dqfunc (void)
{
  DQFuncTestDatum data[] = 
    {{{0.55627,  0.896106, {0.471939, 0.592014, 0.0816613, 0.894271, 0.988155,  0.985155, 0.44964, 0.192222, 0.750061, 0.221445, 0.197233,  0.0397501, 0.497374, 0.707924, 0.220638, 0.930558}, {0.0523232,  0.449555, 0.175595, 0.614461}, {0.734523, 0.205769, 0.619325,  0.718355,  0.262583}, {{3, 0}, {3, 0}, {2, 1}, {3, 1}, {2, 1}, {3, 3}, {1,  3}, {1, 3}},  0.613754}, {{0, 0, 0, 0, -7.5072, -11.7109, -14.5521,  1.72972, -5.0714, 8.27048, 4.9275, 1.74637, 7.9166, 9.40335,  3.4485, 0.606356}, {0, -51.4341, -17.8211, -73.4215}}}, {{0.537663,  0.824084, {0.274428, 0.628599, 0.0880236, 0.631862, 0.524367,  0.407154, 0.890791, 0.592112, 0.0269929, 0.699231, 0.670153,  0.661555, 0.97467, 0.249675, 0.494558, 0.0470933}, {0.240147,  0.0439068, 0.875234, 0.328738}, {0.977564, 0.430153, 0.33757,  0.504654,  0.703136}, {{0, 4}, {1, 3}, {3, 2}, {3, 4}, {3, 3}, {0, 2}, {1,  2}, {3, 4}},  0.801553}, {{-5.98587, -9.24282, -2.16574, -6.84197, -5.34564, -7.45476, -3.36158, 1.36645, 0, 0, 0,  0, -16.8847, -19.2534, -15.6508, -11.049}, {-50.9027, -27.0084,  0, -121.393}}}, {{0.249547,  0.872792, {0.178769, 0.394399, 0.358756, 0.28068, 0.151776,  0.695169, 0.688603, 0.619125, 0.177107, 0.445493, 0.194045,  0.572032, 0.936959, 0.401586, 0.318811, 0.243294}, {0.959396,  0.971434, 0.98124, 0.73864}, {0.25626, 0.16988, 0.731694,  0.865848,  0.0774907}, {{2, 1}, {0, 3}, {3, 3}, {2, 4}, {3, 2}, {1, 1}, {1,  1}, {3, 2}},  0.775481}, {{-1.70401, -4.21377, -8.7955, 0.414968, -3.40429,  13.9662, 3.51599, 0.314584, -3.47507,  4.13685, -2.67809, -13.0238, -6.38273, -10.6948, 0.0717682,  1.13181}, {-15.1463, -12.4059, -50.1646, -25.5531}}}, {{0.372938,  0.585168, {0.925714, 0.0803127, 0.684335, 0.966043, 0.748608,  0.63482, 0.49029, 0.394011, 0.811648, 0.233233, 0.171479,  0.150717, 0.852253, 0.261799, 0.190239, 0.412077}, {0.595993,  0.0919191, 0.458546, 0.546229}, {0.518502, 0.316438, 0.0856078,  0.961061,  0.592788}, {{3, 4}, {1, 2}, {1, 0}, {0, 1}, {3, 1}, {3, 3}, {3,  0}, {3, 4}},  0.236125}, {{-2.48237, 5.28978, 2.93893, 0.337433, 2.20516,  1.13719, 4.11086, 1.4906, 0, 0, 0, 0, -3.88481,  1.99578, -2.91745, -9.55556}, {-12.9988, -16.1749,  0, -93.6054}}}, {{0.401273,  0.995018, {0.84418, 0.601306, 0.910983, 0.601007, 0.0325311,  0.368072, 0.739503, 0.45029, 0.180278, 0.106273, 0.549264,  0.038213, 0.584286, 0.014354, 0.0907185, 0.491984}, {0.0657837,  0.697916, 0.0051107, 0.530923}, {0.472996, 0.461791, 0.603838,  0.535904,  0.628817}, {{1, 3}, {3, 4}, {0, 0}, {3, 2}, {0, 0}, {0, 1}, {3,  0}, {0, 3}},  0.860485}, {{5.85951, 8.14476, 0.226187,  2.86828, -2.61834, -6.14123, -10.6721, 1.92861, 0, 0, 0, 0,  2.28332, -1.00352, 3.13063, -7.95252}, {-51.7999, -26.6651,  0, -84.0703}}}, {{0.692855,  0.934897, {0.596285, 0.492413, 0.953352, 0.484606, 0.416007,  0.38614, 0.404088, 0.446393, 0.831721, 0.371786, 0.31337, 0.95441,  0.765938, 0.67387, 0.308259, 0.423487}, {0.292941, 0.212079,  0.704421, 0.887583}, {0.664125, 0.351593, 0.0115658, 0.952686,  0.0678394}, {{1, 4}, {2, 2}, {2, 2}, {3, 0}, {3, 1}, {1, 4}, {0,  2}, {2, 3}},  0.85918}, {{-2.90385, -3.55401, 0.787199,  0.129673, -4.89043, -6.42421, -8.99361, -12.4993, -14.8481, -15.7623, 1.14265, 0.0977824, 0.807491, 12.6486, 9.86505,  0.54683}, {-7.94176, -58.1002, -53.9005, -31.0039}}}};

  for (unsigned int i = 0; i < sizeof (data) / sizeof (data[0]); ++i)
    {
      float dloglambda[4];
      float dtau[(5-1)*4];

      memset (dloglambda, 0, sizeof (dloglambda));
      memset (dtau, 0, sizeof (dtau));

      dqfunc (data[i].inputs.logalpha,
              5,
              data[i].inputs.tau,
              dtau,
              data[i].inputs.loglambda,
              dloglambda,
              data[i].inputs.pz,
              data[i].inputs.ratings,
              8);

      for (unsigned int k = 0; k < 4; ++k)
        {
          assert (fabs (dloglambda[k] - data[i].desired_output.dloglambda[k]) <=
                  1e-4 * (1 + fabs (dloglambda[k]) + fabs (data[i].desired_output.dloglambda[k])) ||
                  (fprintf (stderr, "%g ?= %g\n", dloglambda[k], data[i].desired_output.dloglambda[k]),
                   0));
        }

      for (unsigned int k = 0; k < (5-1)*4; ++k)
        {
          assert (fabs (dtau[k] - data[i].desired_output.dtau[k]) <=
                  1e-4 * (1 + fabs (dtau[k]) + fabs (data[i].desired_output.dtau[k])) ||
                  (fprintf (stderr, "%g ?= %g\n", dtau[k], data[i].desired_output.dtau[k]),
                   0));
        }
    }
}

int 
main (void) 
{
  srand48 (69);

  test_qfunc ();
  test_dqfunc ();

  return 0;
}
